/****************************************************************
*****************************************************************/

*by Xiaodong Fan, fanxiaodong@gmail.com

#delimit ;

cap log close;
clear all;
drop _all;
set more 1;
pause on;

log using log_MSM_plot9.log, replace;

if (strpos("`c(pwd)'", "_pthealth") > 0) {;
    global gpicsdir "pics_pthealth_sipp";
};
else {;
    global gpicsdir "pics_sipp";
};

global gdtadir ".";

include do_globals.do;

***********************************************;
**  Moments from the sipp data;
** outcome: dta_sipp_Moments.dta;
************************************************;
include do_MSM_infile_sippmoments.do;  // need to include do_globals.do first;

****************************************************************;
**   infile simulated individuals;
**   generate profiles form the simulated individual data;

** outcome: dta_MSM_Indiv.dta, dta_MSM_Profiles.dta;
****************************************************************;
include do_MSM_infile_indiv.do;


/*********************************************************************
    Generate profiles from Simulated individual data;
*********************************************************************/
*define programs;

*plot profiles for each type of het;
cap program drop plot_profiles; 
program define plot_profiles;
   local hettype `1';
   keep if ihet == `hettype';
   local lvxlabel "xlabel(20(10)80,grid)";
   local lvlegr0col1 ring(0) col(1) colgap(*0.75) keygap(*0.5) symxsize(*0.9) symysize(*0.75);

   line lfpr t, lpattern("l" "_") ytitle(${gytitle_lfpr},axis(1)) graphregion(color(white)) ||
     line lnw lnw_fd t, lpattern("-.") yaxis(2) `lvxlabel' legend(pos(${gpos_lfpr}) ring(0) col(1) );
   graph export ${gpicsdir}/eps_type`hettype'_lfpr&lnw.eps,replace;

   line lnw lnw_fd t, lpattern("l" "_") ytitle(${gytitle_lnw},axis(1)) graphregion(color(white)) ||
     line H t, lpattern("-.") yaxis(2) `lvxlabel' legend(pos(7) ring(0) col(1) );
   graph export ${gpicsdir}/eps_type`hettype'_lnw.eps,replace;

   line I I_cond t, lpattern("l" "_") ytitle("Investment") graphregion(color(white)) 
        yscale(range(0 0.15)) ylabel(0(0.03)0.15) `lvxlabel' legend(pos(7) `lvlegr0col1');
   graph export ${gpicsdir}/eps_type`hettype'_I.eps,replace;
   line I I_cond t, lpattern("l" "_") ytitle("Investment",axis(1)) graphregion(color(white)) 
        yscale(range(0 0.15)) ylabel(0(0.03)0.15) ||
     line H t, lpattern("-.") yaxis(2) `lvxlabel' legend(pos(7) `lvlegr0col1');
   graph export ${gpicsdir}/eps_type`hettype'_IH.eps,replace;

   line p50I p50I_cond t, lpattern("l" "_") ytitle("Investment") graphregion(color(white)) 
        yscale(range(0 0.15)) ylabel(0(0.03)0.15) `lvxlabel' legend(pos(7) `lvlegr0col1');
   graph export ${gpicsdir}/eps_type`hettype'_p50I.eps,replace;
   line p50I p50I_cond t, lpattern("l" "_") ytitle("Investment",axis(1)) graphregion(color(white)) 
        yscale(range(0 0.15)) ylabel(0(0.03)0.15) ||
     line p50H t, lpattern("-.") yaxis(2) `lvxlabel' legend(pos(7) `lvlegr0col1');
   graph export ${gpicsdir}/eps_type`hettype'_p50IH.eps,replace;

   if (strpos("`c(pwd)'", "_pt") > 0) {;  // part time;
       line lfpr lfpr_ft lfpr_pt t, lpattern("l" "_" "-.") graphregion(color(white)) 
           `lvxlabel' legend(pos(6) ring(0) col(1) order(1 "LFPR" 2 "FT" 3 "PT"));
       graph export ${gpicsdir}/eps_type`hettype'_lfpr_pt.eps,replace;
       line probV0 probVp t, lpattern("l" "_" "-.") graphregion(color(white)) 
           `lvxlabel' legend(pos(6) ring(0) col(1));
       graph export ${gpicsdir}/eps_type`hettype'_probV0p.eps,replace;
   };

   if (strpos("`c(pwd)'", "health") > 0) {;  //  health;
       foreach jjv in lfpr lnw {;
           line `jjv'1 `jjv'2 `jjv'3 `jjv'4 t, lpattern("l" "_" "-" "-.") graphregion(color(white)) 
               `lvxlabel' legend(pos(7) ring(0) col(1) order(1 "Excellent" 2 "Good" 3 "Bad" 4 "Disabled"));
           graph export ${gpicsdir}/eps_type`hettype'_`jjv'_health.eps,replace;
       };
   };

   if ("`hettype'"=="0") {;
       local lvcond65 if t>=22 & t<=65;
       local lvxlabel65 "xlabel(22 30(10)60 65,grid)";
       line lnw_fd lnw_fd_full t `lvcond65', lpattern("l" "_" "_." "-.") ytitle(${gytitle_lnw}) 
           graphregion(color(white)) `lvxlabel65' legend(pos(6) `lvlegr0col1');
       graph export ${gpicsdir}/eps_type`hettype'_lnw_fd_full_65.eps,replace;

       line lnw lnw_fd lnH_worker lnH t `lvcond65', lpattern("l" "_" "_." "-.") ytitle(${gytitle_lnw}) 
           graphregion(color(white)) `lvxlabel65' legend(pos(6) `lvlegr0col1');
       graph export ${gpicsdir}/eps_type`hettype'_lnwH_65.eps,replace;

       line I I_cond t `lvcond65', lpattern("l" "_") ytitle("Investment") graphregion(color(white)) 
           yscale(range(0 0.15)) ylabel(0(0.03)0.15) `lvxlabel65' legend(pos(7) `lvlegr0col1');
       graph export ${gpicsdir}/eps_type`hettype'_I_65.eps,replace;
       line I I_cond t `lvcond65', lpattern("l" "_") ytitle("Investment",axis(1)) graphregion(color(white)) 
               yscale(range(0 0.15)) ylabel(0(0.03)0.15) ||
           line H t `lvcond65', lpattern("-.") yaxis(2) `lvxlabel65' legend(pos(6) `lvlegr0col1');
       graph export ${gpicsdir}/eps_type`hettype'_IH_65.eps,replace;

       line p50I p50I_cond t `lvcond65', lpattern("l" "_") ytitle("Investment") graphregion(color(white)) 
           yscale(range(0 0.15)) ylabel(0(0.03)0.15) `lvxlabel65' legend(pos(7) `lvlegr0col1');
       graph export ${gpicsdir}/eps_type`hettype'_p50I_65.eps,replace;
       line p50I p50I_cond t `lvcond65', lpattern("l" "_") ytitle("Investment",axis(1)) graphregion(color(white)) 
           yscale(range(0 0.15)) ylabel(0(0.03)0.15) ||
       line p50H t `lvcond65', lpattern("-.") yaxis(2) `lvxlabel65' legend(pos(6) `lvlegr0col1');
       graph export ${gpicsdir}/eps_type`hettype'_p50IH_65.eps,replace;
   };

   line A t if t<${gvNT}, lpattern("l") graphregion(color(white))  ||
      line AIME t, lpattern("-.") yaxis(2) `lvxlabel';
   graph export ${gpicsdir}/eps_type`hettype'_AA.eps,replace;

   line C t if t<${gvNT}, lpattern("l") graphregion(color(white)) `lvxlabel';
   graph export ${gpicsdir}/eps_type`hettype'_C.eps,replace;

   line C C_orig t if t<=${gvNT}, lpattern("l") graphregion(color(white)) `lvxlabel'
        legend(pos(11) ring(0) col(1) order(1 "Adult equiv" 2 "Family"));
   graph export ${gpicsdir}/eps_type`hettype'_C_orig.eps,replace;

   line single mrdspnw mrdspw t if t<=${gvNT}, lpattern("l") graphregion(color(white)) `lvxlabel'
        legend(pos(1) ring(0) col(1) order(1 "Single" 2 "Mrd SP NW" 3 "Mrd SP W"));
   graph export ${gpicsdir}/eps_type`hettype'_FS.eps,replace;

   line sdlnw sdlnw_d t if t<${gvNT}, lpattern("l" "_") graphregion(color(white)) `lvxlabel'
          legend(pos(11) ring(0) col(1) order(1 "Simulation" 2 "Data"));
   graph export ${gpicsdir}/eps_type`hettype'_sdlnw.eps,replace;

   line sdH t if t<${gvNT}, lpattern("l") graphregion(color(white)) `lvxlabel';
   graph export ${gpicsdir}/eps_type`hettype'_sdH.eps,replace;
   line sdA t if t<${gvNT}, lpattern("l") graphregion(color(white)) `lvxlabel';
   graph export ${gpicsdir}/eps_type`hettype'_sdA.eps,replace;
   line sdAIME t if t<${gvNT}, lpattern("l") graphregion(color(white)) `lvxlabel';
   graph export ${gpicsdir}/eps_type`hettype'_sdAIME.eps,replace;

   * transition 1to0 and 0to1;
   foreach iv in lfpr1to0 lfpr0to1  {;
      *local lvtperiod "if t>=${gvdata0} & t<=${gvdata10}";
      if ("`iv'"=="lfpr1to0") {;
         local lvytitle "E to U";
         local lvpos = 11;
      };
      else if ("`iv'"=="lfpr0to1") {;
         local lvytitle "U to E";
         local lvpos = 7;
      };
      line `iv' `iv'_d a`iv' a`iv'_d t if t<${gvNT}, lpattern("l" "_" "l" "_") graphregion(color(white)) `lvxlabel'
            ytitle("`lvytitle'") 
            legend(pos(`lvpos') ring(0) col(1) order(1 "Simulation" 2 "Data" 3 "Sim: Aggregate" 4 "Data: Aggregate"));
      graph export ${gpicsdir}/eps_type`hettype'_`iv'.eps,replace;
   };

end;


/*******************************************************************
  Moments generated in Fortran program on the simulated individuals
********************************************************************/
prog_MSM_moments_infile "txt_MSM_Moments.txt" "";

/****************
* first difference of lnw;
line lnw_f lnw_fd_f t if t<=${gvNT}, lpattern("l" "_" "-" "-.") graphregion(color(white))
       xlabel(${gvdata0} 30(10)60 ${gvdata10},grid)
       ytitle("Simulated wage profiles") 
       legend(pos(6) ring(0) col(1) lab(1 "Log wages") lab(2 "Log wages (FD)"));
graph export ${gpicsdir}/eps_moments_lnw_fd_MSM.eps,replace;
*************/

* merge with Moments from the sipp data;
merge 1:1 t using ${gdtadir}/dta_sipp_Moments.dta;
drop if _merge==2;
drop _merge;
sort t;
save ${gdtadir}/dta_MSM_Moments.dta, replace;

/**********************************************
   Profiles (without _f, changed to _indv) 
   and moments (with _f)
**********************************************/
use dta_MSM_Profiles.dta, clear;
sort t ihet;
merge m:1 t using ${gdtadir}/dta_MSM_Moments.dta, nogenerate;
sort ihet t;
save dta_MSM_ProfilesMoments.dta, replace;
* iv: 0 -- all, 1-9 het type;
forvalues iv = 0/9 {;
   preserve;
      *drop if t>=70;
      foreach jv in lnw {;
          replace `jv' = . if t>=80;
      };    
      plot_profiles `iv';
   restore;
};

* 101--single, 102 married spouse not working, 103--married spouse working;
* forvalues iv = 101/103 {;
*    preserve;
*       plot_profiles `iv';
*    restore;
* };


/*******************************************************************************
  merge with profiles generated from the simulated individual data in STATA
*******************************************************************************/
local lvpms lfpr lfpr_ft lfpr_pt lnw sdlnw lnw_fd lfpr1to0 lfpr0to1 alfpr1to0 alfpr0to1 C  // lnw_fe_full lnw_fe
                lfpr_diff_E2G lfpr_diff_G2B lfpr_diff_B2D lnw_diff_E2G lnw_diff_G2B lnw_diff_B2D
                probV0 probVp ssa irecss lfpr_exit pH41 pH241 plnw41
                p1A p2A p5A p50A p95A p98A p99A p1H p2H p5H p50H p95H p98H p99H
                p50lnw p50inc p50C p50I p50I_cond p50labor;  // p50lnw_fe
*first of all, regular moments;
use ihet t `lvpms' using ${gdtadir}/dta_MSM_Profiles.dta if ihet==0, clear;
drop ihet;
foreach iv in `lvpms' {;
   ren `iv' `iv'_indv;
};
sort t;
merge 1:1 t using ${gdtadir}/dta_MSM_Moments.dta, nogenerate;
label var t "Age";
foreach iv in lfpr lfpr_ft lfpr_pt lnw sdlnw lnw_fd alfpr1to0 alfpr0to1 C plnw41 ssa irecss lfpr_exit  //lnw_fe
              lfpr_diff_E2G lfpr_diff_G2B lfpr_diff_B2D lnw_diff_E2G lnw_diff_G2B lnw_diff_B2D {;
   cap label var `iv'_f "Simulation";
   cap label var `iv'_d "Data";
   cap label var `iv'_indv "Indv_Stata";
};
foreach iv in lfpr1to0 lfpr0to1  {;
   cap label var `iv'_d "Data";
   label var `iv'_indv "Indv_Stata";
};

sort t;
save ${gdtadir}/dta_MSM_Moments.dta, replace;

* plot moments;
local lvmoments0 lfpr lnw sdlnw lnw_fd C ssa plnw41;
if (strpos("`c(pwd)'", "_pt") > 0) {;
    local lvmoments `lvmoments0' lfpr_pt;
    local lvmoments0 `lvmoments';
};
if (strpos("`c(pwd)'", "health") > 0) {;
    local lvmoments `lvmoments0' lfpr_diff_E2G lfpr_diff_G2B lfpr_diff_B2D;
    local lvmoments0 `lvmoments';
};
local lvmoments `lvmoments0';
*local lvmoments lfpr lfpr_pt lnw sdlnw lnw_fd C ssa lfpr_diff_E2G lfpr_diff_G2B lfpr_diff_B2D plnw41;

local lvxlabelNT "xlabel(18 20(10)${gvNT},grid)";

foreach iv of local lvmoments {;
   local lyange "";
   local lylabel "";

   if ("`iv'"=="lfpr_diff_E2G" | "`iv'"=="lfpr_diff_G2B" | "`iv'"=="lfpr_diff_B2D") {;
      local lvtperiod "if t>=${gvdata_diff} & t<=${gvdata10}";
      local lvxlabel "xlabel(${gvdata_diff}(5)60 ${gvdata10},grid)";
   };
   else if ("`iv'"=="ssa") {;
      local lvtperiod "if t>=${gvdata_ss} & t<=${gvdata10_ss}";
      local lvxlabel "xlabel(${gvdata_ss}(1)${gvdata10_ss},grid)";
   };
   else if ("`iv'"=="plnw41") {;
      local lvtperiod "if t>=41 & t<=${gvdata10}";
      local lvxlabel "xlabel(40(5)${gvdata10},grid)";
   };
   else {;
      local lvtperiod "if t>=${gvdata0} & t<=${gvdata10}";
      local lvxlabel "xlabel(${gvdata0} 30(10)60 ${gvdata10},grid)";
   };

   line `iv'_f `iv'_d t `lvtperiod', lpattern("l" "_" "-" "-.") 
         graphregion(color(white)) `lvxlabel'
         ytitle(${gytitle_`iv'}) yscale(range(${gyrange_`iv'})) ylabel(${gylabel_`iv'})
         legend(pos(${gpos_`iv'}) ring(0) col(1) );
   graph export ${gpicsdir}/eps_moments_`iv'.eps,replace;


   * plot for longer period;
   if ("`iv'"=="lfpr_diff_E2G" | "`iv'"=="lfpr_diff_G2B" | "`iv'"=="lfpr_diff_B2D") {;
       local lvtperiod "if t>=${gvdata_diff} & t<=${gvNT}";
       local lvxlabel "xlabel(${gvdata_diff}(5)${gvNT},grid)";
       line `iv'_f `iv'_d t `lvtperiod', lpattern("l" "_" "-" "-.") graphregion(color(white)) `lvxlabel' 
             ytitle(${gytitle_`iv'}) 
             legend(pos(${gpos_`iv'}) ring(0) col(1) );
       graph export ${gpicsdir}/eps_moments80_`iv'.eps,replace;

       local lvtperiod "if t>=${gvdata_diff} & t<=75";
       local lvxlabel "xlabel(${gvdata_diff}(5)70 75,grid)";
       line `iv'_f `iv'_d t `lvtperiod', lpattern("l" "_" "-" "-.") graphregion(color(white)) `lvxlabel' 
             ytitle(${gytitle_`iv'}) 
             legend(pos(${gpos_`iv'}) ring(0) col(1) );
       graph export ${gpicsdir}/eps_moments${gvNT}_`iv'.eps,replace;
   };
   else if inlist("`iv'", "lfpr", "lfpr_pt", "lnw", "lnw_fd", "sdlnw", "C") {;
      local lvtperiod "if t>=18 & t<=${gvNT}";
      local lvxlabel "xlabel(20(10)${gvNT},grid)";
      line `iv'_f `iv'_d t `lvtperiod', lpattern("l" "_" "-" "-.") graphregion(color(white)) `lvxlabel' 
             ytitle(${gytitle_`iv'}) 
             legend(pos(${gpos_`iv'}) ring(0) col(1) );
      graph export ${gpicsdir}/eps_moments80_`iv'.eps,replace;
      local lvtperiod "if t>=${gvdata0} & t<=75";
      local lvxlabel "xlabel(${gvdata0} 30(10)70 75,grid)";
      line `iv'_f `iv'_d t `lvtperiod', lpattern("l" "_" "-" "-.") graphregion(color(white)) `lvxlabel' 
             ytitle(${gytitle_`iv'}) 
             legend(pos(${gpos_`iv'}) ring(0) col(1) );
       graph export ${gpicsdir}/eps_moments${gvNT}_`iv'.eps,replace;
   };
    
   ***** RESULTS: they are identical *********************; 
   line `iv'_f `iv'_indv t if t>=18 & t<=${gvNT}, lpattern("l" "_" "-" "-.") graphregion(color(white))  `lvxlabelNT'
         ytitle(${gytitle_`iv'})   // yscale(range(${gyrange_`iv'})) ylabel(${gylabel_`iv'})
         legend(pos(${gpos_`iv'}) ring(0) col(1) );
   graph export ${gpicsdir}/eps_moments_FortranVsStata_`iv'.eps,replace;
};


line lfpr_exit_indv lfpr_exit_d ssa_indv ssa_d t if t>=50 & t<=80, lpattern("l" "_" "-" "-.") graphregion(color(white))
         xlabel(50(5)80,grid) xline(62 65, lcolor(gs10) lwidth(medthick)) 
         ytitle("Labor Force Exit") legend(pos(1) ring(0) col(1) order(1 "LF exit: Sim" 2 "LF exit: Data" 3 "SSA: Sim" 4 "SSA: Data"));
graph export ${gpicsdir}/eps_moments_lfpr_exit.eps,replace;


foreach iv in lfpr1to0 lfpr0to1  {;
   local lvtperiod "if t>=${gvdata0} & t<=${gvdata10}";
   if ("`iv'"=="lfpr1to0") {;
      local lvytitle "E to U";
      local lvpos = 11;
   };
   else if ("`iv'"=="lfpr0to1") {;
      local lvytitle "U to E";
      local lvpos = 7;
   };

   line a`iv'_f a`iv'_indv t `lvtperiod', lpattern("l" "_" "l" "l") graphregion(color(white))
         xlabel(${gvdata0} 30(10)60 ${gvdata10},grid) ytitle("`lvytitle'") 
         legend(pos(`lvpos') ring(0) col(1) order(1 "Simulation" 2 "Indv_Stata"));
   graph export ${gpicsdir}/eps_moments_FortranVsStata_a`iv'.eps,replace;

   line `iv'_indv `iv'_d a`iv'_f a`iv'_d t `lvtperiod', lpattern("l" "_" "l" "l") graphregion(color(white))
         xlabel(${gvdata0} 30(10)60 ${gvdata10},grid) ytitle("`lvytitle'") 
         legend(pos(`lvpos') ring(0) col(1) order(1 "Simulation" 2 "Data" 3 "Sim: Aggregate" 4 "Data: Aggregate"));
   graph export ${gpicsdir}/eps_moments_`iv'.eps,replace;
};

line lnw_f lnw_d lnw_fd_f lnw_fd_d t if t>=${gvdata0} & t<=${gvdata10}, lpattern("l" "_" "-" "-.") graphregion(color(white))
         xlabel(${gvdata0} 30(10)60 ${gvdata10},grid) ytitle("Log Wages & Log Wages (FD)") 
         legend(pos(6) ring(0) col(1) order(1 "Log Wages: Simulation" 2 "Log Wages: Data" 3 "Log Wages (FD): Simulation" 4 "Log Wages (FD): Data"));
graph export ${gpicsdir}/eps_moments_lnwNFD.eps,replace;

* plot both transitions in one figure;
line alfpr1to0_f alfpr1to0_d alfpr0to1_f alfpr0to1_d t if t>=${gvdata0} & t<=${gvdata10}, lpattern("l" "_" "-" "-.") 
        graphregion(color(white))
        xline(${gvdata0} 62 ${gvdata10}, lcolor(gs10) lwidth(medthick)) xlabel(${gvdata0} 30(10)50 62 ${gvdata10},grid)
        ytitle("LFPR transition") yscale(range(0 0.4)) ylabel(0(0.1)0.4)
        legend(pos(11) ring(0) col(1) rowgap(*0.25) symxsize(*1) order(1 "E to U: Simulation" 2 "E to U: Data" 3 "U to E: Simulation" 4 "U to E: Data"));
graph export ${gpicsdir}/eps_moments_transition.eps,replace;


* model with health, the lfpr_diff graph;
if (strpos("`c(pwd)'", "health") > 0) {;
    line lfpr_diff_E2G_f lfpr_diff_E2G_d lfpr_diff_G2B_f lfpr_diff_G2B_d
         lfpr_diff_B2D_f lfpr_diff_B2D_d t if t>=${gvdata_diff} & t<=${gvdata10}, 
             lpattern("l" "_" "_-" "_.." "-" ".") graphregion(color(white)) xlabel(${gvdata_diff}(5)60 ${gvdata10},grid)
             ytitle(${gytitle_lfpr_diff}) yscale(range(${gyrange_lfpr_diff})) ylabel(${gylabel_lfpr_diff})
             legend(pos(1) ring(0) col(1) order(1 "E-G: simulation" 2 "E-G: data" 
                   3 "G-B: simulation" 4 "G-B: data" 5 "B-D: simulation" 6 "B-D: data") );
    graph export ${gpicsdir}/eps_moments_lfpr_diff.eps,replace;
};

cap log close;



